#include <cmath>
#include <cstdlib>
#include <sstream>
#include <fstream>
#include <string>

#include "sampling.h"

const int R = 32;
const double PI = atan(1.)*4;

void Sampling::norm(const vector<float3>& vtlist, float scale) {
	vector<float3> cv_list(vtlist.size());
	float3 objmax, objmin, cp;
	float avg_dist = 0;

	/* find the "center of mass" of the object */
	objmax = vtlist[0];
	objmin = vtlist[0];

	for (int i = 1; i < vtlist.size(); i++){
		objmax.x = max(objmax.x,vtlist[i].x);
		objmax.y = max(objmax.y,vtlist[i].y);
		objmax.z = max(objmax.z,vtlist[i].z);
		objmin.x = min(objmin.x,vtlist[i].x);
		objmin.y = min(objmin.y,vtlist[i].y);
		objmin.z = min(objmin.z,vtlist[i].z);
	}

	cp = (objmin+objmax) * 0.5;

	/* move the "center of mass" to origin, and
	 * normalize all vertices so that the avg_dist is R/2 */
	for (int i = 0; i < vtlist.size(); i++) {
		cv_list[i] = vtlist[i]-cp;
		avg_dist += sqrt((pow(cv_list[i].x,2) +
						  pow(cv_list[i].y,2) +
						  pow(cv_list[i].z,2)));
	}
	avg_dist /= vtlist.size();

	//ofstream out("test.obj"); /* out put an obj file to see the scaling and transform is correct */
	vertex_list.clear();
	for (int i = 0; i < cv_list.size(); i++) {
		vertex_list.push_back( cv_list[i] * (scale/avg_dist) );
		//out << "v " << nv_list[i].x << " " << nv_list[i].y << " " << nv_list[i].z << " " << endl;
	}
}

void Sampling::objloader(char* filename, int nflag) {
	int polyid;
	string cmd;
	string buffer;
	string pindex;
	float3 tempv;
	polygon tempp;
	vector<polygon> polygon_list;
	vector<float3> vlist_temp;
	ifstream in(filename);
	
	if (!in) {
		cout<<"ERROR: Cannot open file: "<<filename<<endl;
		cout<<"Exit"<<endl;
		exit(1);
	}

	polygon_list.clear();
	cout << "# Reading ";
	cout << filename << " ..." << endl;
	while (getline(in,buffer)) {
		cmd = "#";
		istringstream is(buffer);
		is >> cmd;
		if (cmd=="v") {
			is >> tempv;
			vlist_temp.push_back(tempv);
		}
		else if (cmd=="f") {
			tempp.node.clear();
			while (is>>pindex) {
				size_t slash=0;
				slash=pindex.find("/");
				if (slash>0){
					pindex=pindex.substr(0,slash);
				}
				polyid=atoi(pindex.c_str());
				tempp.node.push_back(polyid);
			}
			polygon_list.push_back(tempp);
		}
		else {
			continue;
		}
	}
	in.close();
	
	int3 tri_temp;
	triangle_list.clear();
	vector<polygon>::iterator itp;
	
	for (itp = polygon_list.begin(); itp<polygon_list.end(); itp++) {
		tri_temp.a = itp->node[0];
		for (int i = 1; i < itp->node.size()-1; i++) {
			tri_temp.b = itp->node[i];
			tri_temp.c = itp->node[i+1];
			triangle_list.push_back(tri_temp);
		}
	}

	if(nflag==0) {
		vertex_list.clear();
		vertex_list=vlist_temp;
	} else if (nflag==1){
		cout<<"## vertices are translated and normalized"<<endl;
		norm(vlist_temp, (R/2));
	}

}

int RayTracer::trace (vector<float3> nv,float theta, float phi, float r) {
	float3 e,d;
	float a3,b3,c3,d3,e3,f3,g3,h3,i3,j3,k3,l3;
	float beta,gama,ist=0,M;

	/* normalized direction "d" */
	d.x = sin(theta)*cos(phi);
	d.y = cos(theta);
	d.z = sin(theta)*sin(phi);
	e.x = 0;
	e.y = 0;
	e.z = 0;
	float3 v[3];

	for (int i = 0; i < triangle_list.size(); i++) {
		v[0] = nv[triangle_list[i].a-1];
		v[1] = nv[triangle_list[i].b-1];
		v[2] = nv[triangle_list[i].c-1];

		a3 = v[0].x - v[1].x;
		b3 = v[0].y - v[1].y;
		c3 = v[0].z - v[1].z;
		d3 = v[0].x - v[2].x;
		e3 = v[0].y - v[2].y;
		f3 = v[0].z - v[2].z;
		g3 = d.x;
		h3 = d.y;
		i3 = d.z;
		j3 = v[0].x-e.x;
		k3 = v[0].y-e.y;
		l3 = v[0].z-e.z;

		M = a3*(e3*i3-h3*f3) + b3*(g3*f3-d3*i3) + c3*(d3*h3-e3*g3);
		beta = (j3*(e3*i3-h3*f3) + k3*(g3*f3-d3*i3) + l3*(d3*h3-e3*g3)) / M;
		gama = (i3*(a3*k3-j3*b3) + h3*(j3*c3-a3*l3) + g3*(b3*l3-k3*c3)) / M;
		ist = (f3*(a3*k3-j3*b3) + e3*(j3*c3-a3*l3) + d3*(b3*l3-k3*c3)) / M;

		if ((beta+gama)<1 && beta>0 && gama>0) {
//			cout << ist << endl;
			if ((-1*ist)>(r-0.5) && (-1*ist)<(r+0.5)) {
				return 1;
			}
		}
	}
	return 0;
}

void RayTracer::getsample(const char* filename, int BW) {
	cout<<"## sampling by ray-tracer"<<endl;
	int tem_sp;
	string temp_fn(filename);
	string::size_type nbegin;
	string::size_type nend;

	nbegin = temp_fn.rfind('/',temp_fn.size());
	nend = temp_fn.rfind('.',temp_fn.size());
	if (nbegin == string::npos || nend == string::npos) {
		cout<<"wrong file name"<<endl;
		exit(0);
	}
	string objname = temp_fn.substr((nbegin+1),(nend-nbegin-1));
	objname = "./data/"+objname+".dat";
	ofstream outf(objname.c_str());
	
	float phi,theta;
	for (float r = 1; r <= R; r++) {
		phi = 0;
		theta = 0;
		outf << "# " << r << endl;
        
		/* samples the surface of radius = r at 2B*2B points */
		for (float k = 0; k < (2*BW); k++) {
			for (float j = 0; j < (2*BW); j++) {
				phi = (PI*k/BW);
				theta = (PI*(2*j+1)/4/BW);
				tem_sp = trace(vertex_list,theta,phi,r);
				outf << tem_sp << endl;
//				outf << "0" << endl;  // required by the origin s2kit of real and imagine parts
			}
		}
	}
	outf.close();
}
